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Abstract 


Direct numerical simulation (DNS) is used to study the mechanism of generation and 
evolution of turbulence structures in a temporally evolving supersonic swirling round jet 
and also to examine the resulting acoustic radiations. Fourier spectral expansions are used 
in the streamwise and azimuthal directions and a 1-D b-spline Galerkin representation 
is used in the radial direction. Spectral-like accuracy is achieved using this numerical 
scheme. Direct numerical simulations, using the b-spline spectral method, are carried 
out starting from mean flow initial conditions which are perturbed by the most unstable 
linear stability eigenfunctions. It is observed that the initial helical instability waves evolve 
into helical vortices which eventually breakdown into smaller scales of turbulence. ‘Rib’ 
structures similar to those seen in incompressible mixing layer flow of Rogers and Moser 1 
are observed. The jet core breakdown stage exhibits increased acoustic radiations. 
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I. Introduction 


The object of this work is to obtain highly accurate solutions of turbulent round jets. 
The solutions obtained help understand the various turbulent scales and mechanisms of 
turbulence generation in the evolution of a compressible round jet. These accurate flow 
solutions will be used to estimate acoustic radiation in the near-field region. There has 
been some work 2-4 in the field of compressible round jets at supersonic Mach numbers 
but none, to the authors’ knowledge, in transition to turbulence in supersonic jets. 

Over the past few decades Direct Numerical Simulations (DNS) have become a very 
powerful tool to obtain highly accurate flow simulations 1,5_9 . In these simulations no 
models are used, the flow equations are solved using highly accurate numerical schemes. 
Hybrid spectral methods form a class of such highly accurate numerical methods which 
have been successfully used by the researchers mentioned above and also by Spalart et al 
10 . For the present work a hybrid b-spline spectral method developed by Moser et al ° is 
used. 

In order to use DNS at high enough Reynolds number to get sufficient turbulent 
structures the temporal jet problem is studied using periodicity in the axial direction. 
This allows for the application of spectral expansion in the axial direction. Physically this 
means that the turbulent structures in the jet are repeated in successive downstream cells 
and grow in time instead of being gradually modified downstream into a jet plume. There 
is an approximate correlation between the growth of structures with time in the temporal 
jet and the growth of structures in the spatially growing jet as one follows the structures 
downstream. Spectral accuracy helps capture smaller turbulent scales at the expense of 
some compromise to the overall jet structure. 

The compressible round jet is simulated by using Fourier expansions in the azimuthal 
and streamwise direction and a 1-D b-spline basis representation in the radial direction. 
The simulation starts with the mean flow, upon which the linear stability eigenfunctions are 
superimposed, and develops temporally. Various stages of growth, evolution and eventual 
breakdown of vortical structures are noted. Preliminary acoustic estimates are presented. 
These give a reasonable picture of the source of most of the turbulent noise and a fair 
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picture of the directivity of these acoustic waves. 


II. Governing Equations 

The compressible Navier-Stokes equations written in cylindrical coordinates in non- 
dimensional form are, 
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p Uk, Re = Pr = Tj is the jet radius, Re , Pr are the 

Reynolds number and Prandtl number, C p is the specific heat at constant pressure, and $ 
is the viscous dissipation given by 
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III. Linear Stability Analysis 
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Inviscid linear stability analysis is performed for jets with varying degrees of swirl for 
a Mach 2 jet. In all the cases studied here it is found that the Kelvin-Helmholtz instability 
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dominates. Helical waves of the form exp{-k z z + rnO) where k z and m are axial and 
azimuthal wavenumbers, are found to be most unstable with m ^ 0. This is similar to 
the compressible shear layer results of Sandham and Reynolds 8,9 , who found that waves 
oblique to the flow direction were most unstable. The main results of the analysis are shown 
in fig 1 which shows growth rate versus wavenumber for a series of swirl ratios. It is seen 
that swirl increases the growth rates of some higher azimuthal disturbances depending on 
their sense of rotation. From the figure it is also evident that positive swirl while increasing 
growth rates also shifts the most unstable waves to higher m and k z , whereas a negative 
swirl shifts the most unstable modes to lower wavenumbers. This is an interesting result 
since the addition of positive swirl would encourage shorter waves to grow faster and on 
the other hand a negative swirl tends to aid the growth of longer waves. 

The results obtained from this analysis are used to compute the most unstable eigen- 
functions, which are used to perturb the mean flow for non-linear DNS calculations. For 
details on the stability analysis refer to Rao 11 . 


IV. Numerical Formulation 


The flow variables are expanded using Fourier series in the two periodic directions, 
viz. the azimuthal (9) and the axial (z) directions. In the non-periodic or radial direction 
(r) 1-dimensional b-splines are used as interpolating functions. 

B-splines 12 of order n are piecewise polynomials of degree n having n — 1 continuous 
derivatives. The n th derivative has a jump at knot points within the interval. Since they 
have a high degree of continuity derivative quantities (like vorticity) can be smoothly and 
accurately represented. B-splines have local support and hence boundary conditions and 
any other conditions which are localized can be easily implemented. Also the local support 
leads to sparse block diagonal matrices which can be efficiently stored and solved. 

The continuity equation, for example, can be written in Galerkin form for any kg and 
k z as 
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where, 

a = ^2<r k (z,6,t) b k (r) , m e = b k (r) , 

k k 

m r = '^T J m r k {z,8,t) b k (r) , m z = ^m Zk (z,0,t) b k {r) 

k k 

The Fourier terms are included in the coefficients of the variables, so cr k (z,9,t) = 
°k(t) E k e E kz exp(ikg9 + ik z z) and so on. 

The derivatives in the 9 and z directions are computed by taking a Fast Fourier 
Transform (FFT) into wave space and multiplying by the appropriate wave numbers (kg 
and k z ), and then an inverse FFT is applied to bring it back to physical space. The 
momentum and energy equations can be written in a similar manner. 

Writing the flow equations as discussed above results in a linear system of coupled 
equations to solve simultaneously at each time step. Since b-splines of order n have local 
support onn + 1 knot (node) intervals this results in a 2n + 1 block banded matrix system, 


M f = R, (5) 

where M is the resulting mass matrix, / is the column vector of nodal values to be solved 
for, and R is a column matrix resulting from the RHS of the governing equations. A low 
storage 3 rd order Runge-Kutta scheme designed by Wray 13 is used for time integration. 

V. Regularity Requirements and Modal Reduction 

In the cylindrical coordinate system the origin (r — 0) is a source of concern since 
some of the functions do not remain analytic. From a mathematical point of view the 
flow variables should be single valued and finite. To enforce this the polynomial expansion 
functions must satisfy some regularity requirements 14 . The z-component of the velocity 
should be represented as, 
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where P 2 (r 2 ; m, A:) is a polynomial in r 2 . Scalars and ^-components of all vectors should 
be represented in a similar manner. The 9 and r-components of the vectors are dependent 
on each other and should be represented as, 

u r (r; m, k) = b(m,k ) P r (r 2 ; m,k) e l m 6 e % kz , 

ug(r; m,k ) = c(m, k) P<?(r 2 ; m,k ) e l m6 e l kz , 

c(m,k ) = ib(m,k ) form > 1, P) 

c(m,k ) = — ib(m,k ) form < — 1, 

P r (0; m, fc) = Pe(0; m,A:) = l, 
for m = 0 6(m, fc) and c(m, A:) are unrelated. 

Another problem arising from a cylindrical mesh is that the aspect ratio of the cells 
near the origin gets very large. The number of azimuthal modes near the origin have to be 
reduced to maintain a good CFL number. This is referred to as mode suppression. The 
number of azimuthal modes are effectively reduced to 2 near the origin and are increased 
successively until the outer boundary has all the azimuthal modes. This reduces the 
accuracy near the axis but increases the allowable time-step, dt, by a significant amount. 

Both of these conditions give rise to a set of constraint equations which are imple- 
mented by replacing some rows in the mass matrix and suitably modifying the RHS vector 
R 7 - n . 


VI. Boundary conditions 

In solving a temporally evolving jet the inflow-outflow boundaries are made periodic 
(see fig 2). This also enables us to use a spectral expansion in the axial direction and 
also takes care of the inflow-outflow boundary conditions. The only boundaries of concern 
are the radial numerical boundaries where a first order non-reflecting boundary condition 


7 



for outflow has been used. This condition prevents incoming waves from infinity by using 
one-dimensional Riemann invariants. This takes the form, 


d , d , 

dt p p™ c °° dt ur 



(8) 


where R is the radius of the outer boundary, and the primed quantities represent perturba- 
tions from the oo quantities, which represent the free stream quantities. The term on the 
right hand side is a correction for a cylindrical boundary. Similar conditions are obtained 
by Engquist and Majda 15,16 and Giles 17 . This boundary condition works well for directly 
incident waves. For the present case the above conditions serve the purpose since the jet 
flow has nearly cylindrical wavefronts incident on the boundary. In addition the outgoing 
wave amplitudes have been mitigated by numerical dissipation due to coarsening of the 
mesh gradually as the outer boundary is approached. 


VII. The DNS of the Swirling Jet 

VII a) Initial Conditions and Computational Domain SetUp 


The simulation is started with a base flow having a half-Gaussian (see fig 3.) axial 
velocity profile with a maximum centerline velocity of Mach 2 and a tangential velocity 
such that there is solid body rotation in the jet core and maximum swirl in the jet shear 
layer diminishing to zero outside of the jet shear layer. This is achieved simply by setting 
Ug = SrU z , where S is a factor used to control the amount of swirl (also called swirl ratio). 
For this simulation S is set to 0.4. The radial mean flow velocity is set to zero. Uniform 
mean temperature is assumed across the jet. The mean pressure P depends on the radial 
coordinate r, and hence so does the mean density p. When the mean flow parameters are 
used to satisfy the Euler equations a relation for P is obtained ; 


dP _ pul 
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Using the above equation in equation 10 it can be shown that 

roo TrUa . 

P(r) = P oc e~Jr rT d \ ( 11 ) 

where T is a constant set to 1, Ug is a known function of r, and P ^ ^ from non- 

dimensionalization. Once P(r) is computed, p can be computed from equation 11. 

The density is non-dimensionalized with p^, and the temperature by Too- Pressure is 
non-dimensionalized by PooO^, where a ^ is speed of sound at infinity. The dimensionless 
mean flow temperature and hence the speed of sound are set to 1. The same mean flow 
conditions are used for the linear stability analysis in section III. The mean flow quantities 
a, fhg, fh z , P are perturbed by the most unstable linear stability eigenfunctions, computed 
in section III (a 1 , u ' e , u' r , u ' z , and p'), giving initial conditions : 

1 , 
a = — + a 

P 

mg —rhg+m'g (= pUg + p'Ug + pu' e ) 
m r — m' r (= pu' r ) 

m z =rh z + m' z (= pU z + p'U z + pu' z ) 

P = P + p' 

The maximum amplitude of these perturbations is taken to be 0.05. 

Though the most unstable eigenvalues for a S = 0.4 occur for m = 5 at k z = 5.5 (fig 
2.), m = 4 and k z = 3 eigenfunctions are used in order to keep the computation similar to 
the non-swirling jet case simulated by Rao 11 . The domain size is chosen to accommodate 
one wavelength of the m = 4 and k z = 3 disturbance wave in both the 6 and z directions, 
so a quarter domain is simulated in the 6 direction and the computational domain length 
in the axial direction is L z = P-. This greatly reduces the size of the simulation and 
allows better resolution for smaller scales. In this case however, the most unstable waves 
have higher frequencies or lower wavelength in the 6 and z directions. Hence due to the 
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fluid physics, more unstable mode waves with higher frequencies can be generated, and if 
they grow faster than the rn = 4 and k z = 3 waves, the domain can exhibit more than 
one wavelength of these shorter wavelength waves. There is thus a possibility for multiple 
structures with this shorter wavelength to evolve, and thus the possibility of capturing 
vortex pairing dynamics if two vortex structures develop in one domain length. This is 
indeed the case as will be seen below in section Vllb. The outer domain is chosen at 4 jet 
radii (R = 4 rj). This is considered far enough not to encounter any mass flow across the 
boundary. The non-reflecting boundary conditions are applied at this boundary. 

The Reynolds number (Re) for the simulation based on jet radius (rj) and speed of 
sound in the core (a) is kept the same at 2500. The simulation is started from a coarse mesh 
which is refined as the simulation progresses. The resolution requirements are controlled 
by monitoring the energy spectra and also by visual inspection of the flow field. It is noted 
that maximum resolution is needed during transition. The final maximum resolution used 
is 128 * 167 * 144 (No *N y * N z ). 128 Fourier modes represent a quadrant in the azimuthal 
direction, 167 b-splines in the radial direction from 0 - R, and 144 Fourier modes in the 
axial direction over L z are used. The mesh is non-uniform in the radial direction with 
more clustering in the jet shear layer. 

VII b) The DNS Results 

The evolution of the swirling jet is shown in figures 4a and 4b as a reconstruction of the 
temporal jet, constructed by placing temporal sections from different times adjacent to each 
other. The convective velocity of the large scale structures is used for the reconstruction. 
The density plot (fig. 4a) clearly outlines the different stages of evolution. The initial 
frames show a uniformly varying density from the axis to the shear layer, until t = 5.2 
which marks the commencement of core breakdown. Also clearly visible is the edge of the 
jet shear layer. The last two frames show flow without much of a core. Some Mach waves 
can be seen in the last few frames after core breakdown (seen more clearly in a similar 
pressure plot shown below in fig 14). 

The vorticity magnitude plot (fig. 4b) illustrates an initially sharp shear layer with 
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helical instability waves developing into helical vortices. There is distinct formation of 
roller at time t = 2.6. Since the domain size is such that it supports a wavelength of the 
disturbance smaller than its length (explained above in Vila) there is the formation of two 
rollers (also seen in the r - 9 cross-section plot fig. 9a). These rollers eventually pair and 
this pairing mechanism can be clearly seen at t = 3.9. The pairing is complete by t = 5.2 
leaving only one large scale structure, and by t = 6.5 this structure rolls up into a large 
roller. At t = 7.8 the roller has broken into smaller more intense vortices. 

The energy spectrum is plotted versus the axial (fig. 5a) and the azimuthal (fig. 5b) 
wave numbers at different times. The solid line shown is the A; -5 / 3 spectrum. A broad 
energy spectrum is observed even though the initial conditions had energy only in the k z 
= 1 and kg — 1 modes. This indicates that the energy has been cascaded to other scales 
as the jet developed. 

Figure 6, shows variation of mean axial velocity with radius over time. This is used 
to monitor the growth in the mean thickness of the shear layer at different times of the 
simulation. The shear layer which was initially about ().2rj reaches 2. Or 3 at core collapse. 

Since only a m = 4 mode is used to perturb the mean flow, the ensuing flow will have 
a quarter domain symmetry in the azimuthal direction. Thus flow is computed only in 
one quadrant. The other quadrants would have identical flow fields. Hence we plot only 
one quadrant. The instability waves are tracked from initial time through a sequence of 
vorticity iso-surface plots with a few cross-sectional plots for emphasis. 

At t = 0 (fig 7), the iso-surface plots show a flat vortex sheet with a strong helical 
wave instability on it. By t = 1.3 (fig 8) the development of multiple helical disturbances 
is clearly visible. This indicates the start of the evolution of the shorter wavelength dis- 
turbances. This is when the faster growing instability waves (higher frequency) begin to 
dominate the flow. The higher vorticity (u> = 12) is tube-like and sheathed between the 
lower vorticity (w = 6) sheets. These waves become more pronounced at t = 2.6 (fig 9b). 
There is a distinct roll-up of the helical instability waves into helical vortex tubes. Two 
vortex rollers both in the r - 9 (fig 9a) and r-z cross-sections ( see fig 4b) are clearly seen. 
This allows the phenomenon of vortex-pairing. The vortex sheets in the region between the 
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two rollers (this has been called the braid region by some researchers) exhibit secondary 
instabilities. Pairing can be seen between the two helical vortex rollers at t = 3.9 (fig 10a 
- 10b and 4b). These rollers begin merging into one large scale structure. The vorticity is 
still concentrated in long vortex-sheet like structures. Intense vortex tubes or ‘ribs’, which 
are oriented almost perpendicular(locally) to and wrap around the primary helical vortex, 
can initiate in the braid region between the two rollers. (Rogers and Moser 6 , noted very 
intense round ‘ribs’ or streamwise vortices wrapping around large scale spanwise rollers in 
their incompressible mixing layer flow.) These ‘ribs’ develop circular cross-sections as they 
intensify. 

By t = 5.2 (fig lla-llc) the pairing mechanism is complete and the large scale structure 
has started breaking down into smaller more intense vortical structures. This is a rapid 
transition region and is characterized by large acoustic radiations, as will be shown in the 
next section. The r — 6 cross-section shows formation of spiral S-shaped vortex cores along 
the outer edge of the shear layer. These intense vortex cores roll the fiat vortex sheets 
around them. These are ‘rib’ cross-section signatures. The ‘ribs’ are clearly visible in the 
iso-surface plots (u> — 6) as the two large vortex tubes wrapping around the primary vortex 
tube in the opposite sense of the helix f. The primary helical vortex is seen more clearly 
in the u> = 12 plot (fig 11c). 

No more vortex sheet like structures are visible at t = 6.5 (fig 12a-c). Prom fig. 4b, 
it is noted that this was the time of roll-up into a large structure. Most of the vortical 
structures are large, tubular and oriented randomly. The ribs are still persistent, especially 
the one in the center of the plot, it can be seen at all three vorticity levels, indicating that 
the ‘ribs’ have become more intense after getting stretched. At t = 7.8 (fig. 13a-d), large 
scale structures are no longer seen; most of the vorticity is concentrated in smaller more 
intense tubes which are randomly oriented. 

VIII. Acoustic Estimates 

Fig 14 shows the various stages of jet evolution in a reconstruction similar to fig 4. 

f The primary vortex goes from top left to bottom right of the figure and loops around 

the jet core such that the axis of the helix is in the streamwise direction 


12 



Prom this pressure plot it can be seen that the uniformly varying jet core begins to collapse 
at t = 5.2, after the vortex pairing has occurred. The pairing mechanism seems to induce 
early transition. The core breakdown is seen to be the primary source of noise. In the next 
couple of frames Mach waves, which originate at core breakdown, can be seen propagating 
outwards. 

Plotted in fig 15 is the total acoustic energy flux per unit area ( p ' * u r ) radiated 
through a quarter domain of the cylindrical shell at r = 3. This gives a rough estimate of 
when the acoustic energy is radiated. From this plot one sees a sharp rise in the energy 
flux near t = 2. This is when the initial disturbances in the shear layer propagate to 
r = 3. Also around t = 3 marks the growth and non-linear interactions of instability 
waves. The energy flux increases with time indicating that acoustic emissions increase as 
the helical waves grow in amplitude and interact with each other. At t = 7.2 a peak is 
clearly noticed, this corresponds to the fluid physics in the core around t = 4.2 - 5. So as 
expected the vortex pairing and the rapid breakdown of the large scale structures at core 
collapse cause a rise in acoustic radiations. This is similar to results seen by Mitchell et al 
18 , and is consistent with the premise that rapidly varying vortical structures cause most 
of the noise. The ensuing small scale turbulence should give higher frequency noise. At 
core breakdown Mach wave radiations cause most of the acoustic radiation, and this will 
be captured by the pressure energy flux plot around t = 8.2. 

The fluctuation of p over time at a point on the r = 3 shell is plotted in fig 16. It can 
be seen that the frequency of the oscillations increases over time as the flow becomes more 
turbulent. 

IX. Conclusions 

The instability waves are tracked in time as they evolve. Large scale vortex pairing 
is observed which appears to hasten transition to turbulence. Mach waves are observed 
around the time of core collapse. The energy spectra indicate that the energy which was 
initially in one mode has cascaded to higher wave numbers or smaller structures. ‘Rib’ 
structures similar to those seen in incompressible mixing layer flow of Rogers and Moser 
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1 are observed. This is interesting since the flow initial and boundary conditions and the 
governing flow equations of these flows is quite different from the mixing layer flow. 

The acoustic emissions for the jet flow are estimated. It is seen that the core collapse 
of the jet causes increased acoustic emissions. The time series plot of the pressure pertur- 
bation at r = 3 shows most of the fluctuations to be low frequency waves. Slightly higher 
frequency is seen as turbulence develops; smaller vortices give higher frequency emissions. 
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Fig. 1 : Effect of swirl on linear stability growth rates 

Fig. 2 : Schematic of the computational domain 

Fig. 3 : half-Gaussian velocity profile 

Fig. 4a : The reconstructed swirling jet : density 

Fig. 4b : The reconstructed swirling jet : vorticity magnitude 

Fig. 5a : Plots showing the energy spectrum vs axial wavenumber 

Fig. 5b : Plots showing the energy spectrum vs azimuthal wavenumber 

Fig. 6 : Plot showing the variation of axial velocity over time 

Fig. 7 : Total vorticity iso-surface plot at time t = 0 

Fig. 8 : Total vorticity iso-surface plot at time t = 1.3 

Fig. 9a : Total vorticity r — 6 plot at time t = 2.6 

Fig. 9b : Total vorticity iso-surface plot at time t = 2.6 

Fig. 10a : Total vorticity r — 6 plot at time t = 3.9 

Fig. 10b : Total vorticity iso-surface plot at time t = 3.9 

Fig. 11a : Total vorticity r — 6 plot at time t = 5.2 

Fig. lib : Total vorticity (u> = 6) iso-surface plot at time t = 5.2 

Fig. 11c : Total vorticity ( ui = 12) iso-surface plot at time t = 5.2 

Fig. 12a : Total vorticity (oj = 6) iso-surface plot at time t = 6.5 

Fig. 12b : Total vorticity (u = 12) iso-surface plot at time t = 6.5 

Fig. 12c : Total vorticity (u = 16.5) iso-surface plot at time t = 6.5 

Fig. 13a : Total vorticity (u) = 6) iso-surface plot at time t = 7.8 

Fig. 13b : Total vorticity (a; = 12) iso-surface plot at time t = 7.8 

Fig. 13c : Total vorticity (u> = 16.5) iso-surface plot at time t = 7.8 

Fig. 13d : Total vorticity (u> = 20) iso-surface plot at time t = 7.8 

Fig. 14 : The evolution of pressure over time 

Fig. 15 : The acoustic energy radiated over time for the swirling jet 
Fig. 16 : Pressure variation with time for the swirling jet 
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(a) Energy spectrum va k. 
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(b) Energy spectrum vs k e 
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Total vorticlty magnitude ( t = 2.6) 
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25.9982 
23.1095 
20.2208 
17.3321 
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Total vorticlty magnitude ( t= 3.9) 
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2.97346 


Fig. 10a 



Fig. 10b 


Ram M Rao (Physics of Fluids) 





Total vortldty magnitude ( t = 5.2 ) 



1 43.4232 
40.5284 
37.6335 
34.7386 
31.8437 
28.9488 
26.0539 
23.1591 
20.2642 
17.3693 
14.4744 
11.5795 
8.68465 
5.78977 
2.89488 


Fig. 11a 



Fig. lib 


Ram M Rao (Physics of Fluids) 








t * 6.5 



Fig. 12a 


Ram M Rao (Physics of Fluids) 




Fig. 12b 



Ram M Rao (Physics of Fluids) 






I * 7.8 



Fig. 13a 


Ram M Rao (Physics of Fluids) 





Fig. 13b 


Ram M Rao (Physics of Fluids) 





